%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Figure 2 in the paper "Optimal 
%%% Portfolio Choice with Fat Tails and Parameter Uncertainty",
%%% Journal of Financial and Quantitative Analysis (2024),
%%% Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Exponential integral and risk-aversion coefficient
Expfun = @(n,x) integral(@(t)exp(-x.*t)./(t.^n),1,inf);
gam = 1;

%%% Initialize vectors
nuvec = 3:0.1:20;
EUrho01theta03 = zeros(length(nuvec),1);
EUrho03theta03 = zeros(length(nuvec),1);
EUrho05theta03 = zeros(length(nuvec),1);
EUrho01theta05 = zeros(length(nuvec),1);
EUrho03theta05 = zeros(length(nuvec),1);
EUrho05theta05 = zeros(length(nuvec),1);

%%% \theta = 0.3
theta2 = 0.3^2;
Umaxtheta03 = theta2/(2*gam);
j = 1;
for nu = nuvec
    %\rho = 0.1
    rho = 0.1;
    y = @(x) (nu-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nu/(nu-2)));
    varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
    EUrho01theta03(j,1) = (2*eta/(1-rho)-varphi/(1-rho)^3)*Umaxtheta03...
                          -eta*rho/(2*gam*(1-rho)^3);
    %\rho = 0.3
    rho = 0.3;
    y = @(x) (nu-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nu/(nu-2)));
    varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
    EUrho03theta03(j,1) = (2*eta/(1-rho)-varphi/(1-rho)^3)*Umaxtheta03...
                          -eta*rho/(2*gam*(1-rho)^3);
    %\rho = 0.5
    rho = 0.5;
    y = @(x) (nu-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nu/(nu-2)));
    varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
    EUrho05theta03(j,1) = (2*eta/(1-rho)-varphi/(1-rho)^3)*Umaxtheta03...
                          -eta*rho/(2*gam*(1-rho)^3);
    j = j+1;
end

%%% \theta = 0.5
theta2 = 0.5^2;
Umaxtheta05 = theta2/(2*gam);
j = 1;
for nu = nuvec
    %rho = 0.1
    rho = 0.1;
    y = @(x) (nu-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nu/(nu-2)));
    varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
    EUrho01theta05(j,1) = (2*eta/(1-rho)-varphi/(1-rho)^3)*Umaxtheta05...
                          -eta*rho/(2*gam*(1-rho)^3);
    %rho = 0.3
    rho = 0.3;
    y = @(x) (nu-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nu/(nu-2)));
    varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
    EUrho03theta05(j,1) = (2*eta/(1-rho)-varphi/(1-rho)^3)*Umaxtheta05...
                          -eta*rho/(2*gam*(1-rho)^3);
    %rho = 0.5
    rho = 0.5;
    y = @(x) (nu-2).*rho.*x./(2*(1-rho));
    fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
    eta = fzero(fun,0.5*(1+nu/(nu-2)));
    varphi = 2*eta^2*(1-rho)/(nu-eta*(nu-2));
    EUrho05theta05(j,1) = (2*eta/(1-rho)-varphi/(1-rho)^3)*Umaxtheta05...
                          -eta*rho/(2*gam*(1-rho)^3);
    j = j+1;
end

%%% Create Figure 2
fig = figure;
%Subplot 1
subplot(2,2,1)
plot(nuvec,Umaxtheta03.*ones(length(nuvec),1),'Color',[0.47 0.67 0.19],'LineWidth',2);
hold on
plot(nuvec,EUrho01theta03,'-.b','LineWidth',2);
plot(nuvec,EUrho03theta03,'--r','LineWidth',2);
plot(nuvec,EUrho05theta03,':k','LineWidth',2);
hold off
legend('$U(w^\star)$','$U(\hat w)$ $(\rho=0.1)$','$U(\hat w)$ $(\rho=0.3)$',...
       '$U(\hat w)$ $(\rho=0.5)$','Location','southeast','interpreter','latex');
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\theta=0.3$','interpreter','latex');
xlabel('$\nu$','interpreter','latex');
xlim([3 20]);
ylim([-6 0.3]);
xticks([3 5 10 15 20]);
yticks([-6 -4 -2 0]);
%Subplot 2
subplot(2,2,2);
plot(nuvec,Umaxtheta05.*ones(length(nuvec),1),'Color',[0.47 0.67 0.19],'LineWidth',2);
hold on
plot(nuvec,EUrho01theta03,'-.b','LineWidth',2);
plot(nuvec,EUrho03theta03,'--r','LineWidth',2);
plot(nuvec,EUrho05theta03,':k','LineWidth',2);
hold off
set(gca,'TickLabelInterpreter','latex');
grid on
title('$\theta=0.5$','interpreter','latex');
xlabel('$\nu$','interpreter','latex');
xlim([3 20]);
ylim([-10 0.5]);
xticks([3 5 10 15 20]);
yticks([-10 -5 0]);
fontsize(fig,16,"points");

